      PROGRAM TRAJEC(INPUT,OUTPUT,TAPE5=INPUT)
C THIS IS THE DRIVER FOR THE TRAJEC ELECTRON TRAJECTORY CODE. A DESCRIP-
C TION OF THE CODE MAY BE HAD BY RUNNING THE CODE IN ROOKIE MODE, OR
C EXAMINING THE FORMAT STATEMENTS AT THE END OF THIS ROUTINE.
C T.C. STEPHENS > T.N. DELMER, SCIENCE APPLICATIONS INC., MARCH, 1978
C
      DIMENSION SCORE(10,10),XIMAG(100),YIMAG(100),QIMAG(100),XINPT(5)
      COMMON/IMAG/X(3),U(3),T,NSTEP,XERR,YERR,ZERR
      LOGICAL ROOKIE
      INTEGER YES,LINE,POINT,EXIT,ANS
      DATA YES,LINE,POINT,EXIT/"Y","L","P","E"/
C
      PRINT 100
      CALL GETCRD(ANS)
      ROOKIE=(ANS.EQ.YES)
      IF(.NOT.ROOKIE) GO TO 1
      PRINT 300
      CALL GETCRD(ANS)
      IF(ANS.NE.YES) GO TO 1
      PRINT 400
      PRINT 401
1     PRINT 500
      CALL GETCRD(ANS)
      IF(ANS.EQ.EXIT) STOP
      IF(ANS.EQ.LINE) GO TO 2
      IF(ANS.EQ.POINT) GO TO 5
C
C RUN ONE PARTICLE AT A TIME (OPTIONALLY CHANGE INTERNAL PARAMETERS)
C
      CALL PARAM(ROOKIE)
14    PRINT 1300
      IF(ROOKIE) PRINT 600
      CALL GET(6,X)
      CALL IMAGE
      PRINT 1400,X,XERR,YERR,ZERR,T,NSTEP
      CALL GETCRD(ANS)
      IF(ANS.EQ.YES) GO TO 14
      GO TO 1
C
C GENERATE A LINE OF PARTICLES
C
2     PRINT 1500
      IF(ROOKIE) PRINT 700
      CALL GET(5,XINPT)
      NPART=XINPT(1)
      X1=XINPT(2)
      Y1=XINPT(3)
      X2=XINPT(4)
      Y2=XINPT(5)
      DY=(Y2-Y1)/(NPART-1)
      DX=(X2-X1)/(NPART-1)
      DO 4 I=1,NPART
        XX=X1+(I-1)*DX
        YY=Y1+(I-1)*DY
        DO 3 J=1,3
          X(J)=0.
3         U(J)=0.
        X(1)=XX
        X(2)=YY
        CALL IMAGE
        PRINT 1600,XX,YY,X(1),X(2)
4       CONTINUE
      GO TO 1
C
C SOURCE POINT OF MONTE CARLO SELECTED PARTICLES
C
5     PRINT 1700
      IF(ROOKIE) PRINT 800
      CALL GET(4,XINPT)
      NPART=XINPT(1)
      X1=XINPT(2)
      Y1=XINPT(3)
      Z1=XINPT(4)
      IF(NPART.GT.100) GO TO 5
      CALL SOURCE(0,Q,U(1),U(2),U(3))
      DO 7 I=1,NPART
        X(1)=X1
        X(2)=Y1
        X(3)=Z1
        CALL SOURCE(1,Q,U(1),U(2),U(3))
        CALL IMAGE
        XIMAG(I)=X(1)
        YIMAG(I)=X(2)
        QIMAG(I)=Q
        IF(I.NE.1) GO TO 6
        XMIN=X(1)
        XMAX=X(1)
        YMAX=X(2)
        YMIN=X(2)
6       XMIN=AMIN1(XMIN,X(1))
        XMAX=AMAX1(XMAX,X(1))
        YMAX=AMAX1(YMAX,X(2))
        YMIN=AMIN1(YMIN,X(2))
7       CONTINUE
      DO 8 I=1,10
        DO 8 J=1,10
8         SCORE(I,J)=0.
      DO 13 I=1,NPART
        DO 9 J=1,10
          IXBIN=J
          IF(XIMAG(I).LE.(XMIN+J*(XMAX-XMIN)/10.)) GO TO 10
9         CONTINUE
10      DO 11 J=1,10
          IYBIN=J
          IF(YIMAG(I).LE.(YMIN+J*(YMAX-YMIN)/10.)) GO TO 12
11        CONTINUE
12      SCORE(IXBIN,IYBIN)=SCORE(IXBIN,IYBIN)+QIMAG(I)
13      CONTINUE
      PRINT 1800,XMIN,XMAX,YMIN,YMAX,SCORE
      GO TO 1
C
100   FORMAT(" ARE YOU A ROOKIE <<")
300   FORMAT(" WOULD YOU LIKE AN INTRODUCTION <")
400   FORMAT(" THIS IS THE _TRAJEC_ CODE...IT IS DESIGNED TO CALCULATE
     *"/" ACCURATE TRAJECTORIES OF ELECTRONS IN STATIC 3-D ELECTRO-
     *"/" MAGNETIC FIELDS. THE FIELDS ARE UNIFORM BY DEFAULT, BUT MAY
     *"/" BE OVERRIDDEN BY A USERS OWN SET OF FIELDS THROUGH A SIMPLE
     *"/" INTERFACE. GEOMETRY IS CARTESIAN > ALL UNITS ARE CGS (YES,
     *"/" STATVOLTS!). TRAJEC ASSUMES THAT EACH ELECTRON EVENTUALLY
     *"/" ENCOUNTERS AN X,Y _IMAGE PLANE_ AT WHICH TIME SCORING OR OTHER
     *"/" END CLASSIFICATION TAKES PLACE. TRAJEC OPERATES IN 3 MODES!
     *"/"     1)SINGLE PARTICLE, WHERE ONE ELECTRON IS LAUNCHED WITH AN
     *"/"       ARBITRARY INITIAL POSITION > VELOCITY. THIS MODE PERMITS
     *"/"       USER ACCESS TO INTERNAL CONVERGENCE, GEOMETRY > CONTROL
     *"/"       PARAMETERS AS WELL AS EM CHARACTERISTICS.
     *"/"     2)LINE OF PARTICLES, WHERE A LINEAR ARRAY OF ELECTRONS
     *"/"       (ZERO INITIAL VELOCITIES) ACCELERATE TO A FINAL CONFIG-
     *"/"       URATION WHICH IS LISTED.")
401   FORMAT("     3)POINT SOURCE, WHERE MONTE CARLO SAMPLING OF USER
     *"/"       DEFINED ENERGY > ANGLE DISTRIBUTIONS (SUBROUTINE SOURCE)
     *"/"       PRODUCES ELECTRONS FROM SOME POINT WHICH FORM A DIFFUSE
     *"/"       AREA ON THE IMAGE PLANE.
     *"//" IT IS RECOMMENDED THAT THE SINGLE PARTICLE CASE BE RUN FIRST
     *"/" TO ENSURE THAT THE PROBLEM IS CORRECTLY INITIALIZED.
     *"//" AS MENTIONED ABOVE, A USER MAY INTRODUCE HIS OWN TAILOR MADE
     *"/" ELECTRIC POTENTIAL SIMPLY BY REPLACING SUBROUTINE NUPOT(MODE,
     *"/" X,PHI) WITH HIS OWN VERSION. THIS ROUTINE SUPPLIES THE ELEC-
     *"/" TRIC POTENTIAL, PHI, AT POINT X(I) TO THE TRAJEC CODE WHEN
     *"/" CALLED WITH MODE=1. IT IS INITIALIZED (MODE=0) DURING SINGLE
     *"/" ELECTRON RUNS WHERE THE NON UNIFORM POTENTIALS SWITCH HAS BEEN
     *"/" THROWN. AT THAT TIME, THE USER MAY REQUEST SPECIAL INPUTS OR
     *"/" DO SOME SET UP WORK FOR HIS OWN FIELDS. SIMILARLY, NON-UNIFORM
     *"/" MAGNETIC FIELDS ARE SET UP IN SUBROUTINE NUBFLD(MODE,X,B)")
500   FORMAT(//" WOULD YOU LIKE TO RUN A SINGLE PARTICLE, A POINT SOUR
     *CE, A LINE SOURCE, OR EXIT<")
600   FORMAT(" EXAMPLE! 0, 1, 3.5,    3.1E7, 0, 0")
700   FORMAT(" EXAMPLE! 5, -1.3, 1, 1.3, 1")
800   FORMAT(" EXAMPLE! 50, -1, 2.7, 0")
1300  FORMAT(" INPUT STARTING POSITION AND VELOCITY COMPONENTS")
1400  FORMAT(" THE FINAL PARTICLE POSITION IS (X,Y,Z)="3F10.4
     */" ACCUMULATED (WORST CASE) ERRORS IN X,Y,Z ="3F7.2," MICRONS"
     */" TRAVERSE TIME="E9.3," SECONDS OVER "I3," STEPS"
     */" WOULD YOU LIKE TO RUN ANOTHER ONE IN THESE FIELDS<")
1500  FORMAT(" INPUT THE # OF PARTICLES AND THE ENDPOINTS OF THE DESIR
     *ED LINE AS X1,Y1,X2,Y2")
1600  FORMAT(" SOURCE X,Y="2F8.4,"    IMAGE X,Y="2F8.4)
1700  FORMAT(" INPUT # PARTICLES (CURRENT LIMIT 100) AND SOURCE POINT X
     *,Y,Z")
1800  FORMAT(" THE DIFFUSE IMAGE IS REPRESENTED BY A 10X10 GRID OF ACC
     *UMULATED WEIGHTS"//" THE GRID LIMITS ARE X="2F8.4," Y="2F8.4/,
     *" THE GRID IS (X INCREASING TO RIGHT, Y INCREASING DOWN)!"
     *//(10F6.2))
      END
      SUBROUTINE PARAM(ROOKIE)
C
C THIS ROUTINE DEFINES DEFAULT PARAMETER VALUES > ACCEPTS CHANGES.
C T.C. STEPHENS
C
      COMMON/GUTS/DIAG(5),B(3),V,VUNIF,BUNIF,DELTXI,RLIM,ZLIM,HOPMAX,
     * EP,DIMAG,STPTIM,DTIME
      LOGICAL VUNIF,DIAG,BUNIF
      LOGICAL ROOKIE
      INTEGER YES,ANS
      DIMENSION MATCH(66),IFLD(7),UNI(2),TRK(2)
      DATA ZLIM,RLIM,HOPMAX,EP,DELTXI,DIMAG/15.,1.75,3.,.1,.05,.0001/
      DATA STPTIM,DTIME/1000.,.23/,B,V,VUNIF,BUNIF,DIAG,YES/
     * 5.7785,1.8455,112.1,83.333,2*.TRUE.,5*.FALSE.,"Y"/
      DATA MATCH/16,2,"B","(",1,"V",2,"D","E",2,"R","L",2,"Z","L",
     * 2,"H","O",2,"E","P",2,"D","I",2,"S","T",2,"D","T",7,"U","N",
     * "I","F","O","R","M",10,"N","O","N","U","N","I","F","O","R","M"
     * ,7,"T","R","A","C","K","O","F",2,"T","R",2,"L","I",2,"G","O"/
      DATA IFLD/2,2,"P","O",2,"B","F"/
      DATA UNI/4H NON,4H    /,TRK/" ",3HOFF/,IV,IB,ITK/2,2,2/
C
      PRINT 600,ZLIM,RLIM,HOPMAX,EP,DELTXI,DIMAG,B,V,STPTIM,DTIME,
     * UNI(IV),UNI(IB),TRK(ITK)
      IF(ROOKIE) PRINT 700
      IF(ROOKIE) CALL GETCRD(ANS)
      IF(ANS.EQ.YES.AND.ROOKIE) PRINT 800
      IF(ANS.EQ.YES.AND.ROOKIE) PRINT 801
      PRINT 900
1     CALL GETCRD(DUM)
      CALL SCAN(MATCH,IT)
      IF(IT.EQ.0) PRINT 200
      IF(IT.EQ.0) GO TO 1
      IF(IT.EQ.16) GO TO 3
      IF(IT.GT.10) GO TO 2
      IF(IT.EQ.1) IDX=IBCD(DUM)
      CALL SCANC("=")
      X=FPBCD(DUM)
      IF(IT.EQ.1) B(IDX)=X
      IF(IT.EQ.2) V=X
      IF(IT.EQ.3) DELTXI=X
      IF(IT.EQ.4) RLIM=X
      IF(IT.EQ.5) ZLIM=X
      IF(IT.EQ.6) HOPMAX=X
      IF(IT.EQ.7) EP=X
      IF(IT.EQ.8) DIMAG=X
      IF(IT.EQ.9) STPTIM=X
      IF(IT.EQ.10) DTIME=X
      GO TO 1
2     IF(IT.LT.13) CALL SCAN(IFLD,IF)
      IF(IT.EQ.11.AND.IF.EQ.1) VUNIF=.TRUE.
      IF(IT.EQ.11.AND.IF.EQ.2) BUNIF=.TRUE.
      IF(IT.EQ.12.AND.IF.EQ.2) BUNIF=.FALSE.
      IF(IT.EQ.12.AND.IF.EQ.1) VUNIF=.FALSE.
      IF(IT.EQ.13) DIAG(1)=.FALSE.
      IF(IT.EQ.14) DIAG(1)=.TRUE.
      IV=1
      IB=1
      ITK=2
      IF(VUNIF) IV=2
      IF(BUNIF) IB=2
      IF(DIAG(1)) ITK=1
      IF(IT.EQ.15) PRINT 600,ZLIM,RLIM,HOPMAX,EP,DELTXI,DIMAG,B,V,
     * STPTIM,DTIME,UNI(IV),UNI(IB),TRK(ITK)
      GO TO 1
3     IF(DELTXI.GT.(2./5.*RLIM)) PRINT 1000
      IF(.NOT.VUNIF) CALL NUPOT(0,X,PHI)
      IF(.NOT.BUNIF) CALL NUBFLD(0,X,B)
      RETURN
C
200   FORMAT(" UNRECOGNIZED INPUT...TRY AGAIN")
700   FORMAT(" THESE VARIABLES CONTROL INTERNAL CONVERGENCE, GEOMETRY,
     *"/" > UNIFORM E>M POTENTIALS. DO YOU NEED HELP WITH THEM <")
600   FORMAT(/" THE PROGRAM PARAMETERS ARE!"/" ZLIM="G13.6," CM"/" RLIM=
     *"G13.6," CM"/" HOPMAX="G13.6," CM"/" EP="G13.6," MICRONS"/" DELTXI
     *="G13.6," CM"/" DIMAG="G13.6," CM"/" B(I)="3G13.6," GAUSS"/" V="G1
     *3.6," STATVOLT"/" STPTIM="G13.6," NS"/" DTIME="G13.6," NS"/A4,"UNI
     *FORM POTENTIAL"/A4,"UNIFORM B FIELD"/" TRACK "A3)
800   FORMAT(" ZLIM--POSITION ON Z AXIS OF X-Y _IMAGE PLANE_
     *"/" RLIM--RADIAL BOUNDARY OF E>M FIELDS.
     *"/" HOPMAX--MAXIMUM DISTANCE THAT THE ELECTRON IS ALLOWED TO
     *"/"  MOVE PER TIME STEP (TO INSURE THAT EM FIELD DETAIL IS _SEEN_
     *"/"  BY THE PARTICLE).
     *"/" EP--ACCEPTABLE ERROR OR TOLERANCE IN THE PARTICLE POSITION
     *"/"  OVER EACH TIME STEP (IN MICRONS). THE ACCUMULATED ERROR
     *"/"  OVER THE WHOLE TRAJECTORY IS PRINTED AT THE END OF EACH
     *"/"  SINGLE PARTICLE RUN.
     *"/" DELTXI--COORDINATE STEPSIZE IN POTENTIAL DERIVATIVE CALC.
     *"/" DIMAG--Z DISTANCE FROM IMAGE PLANE WHERE ELECTRON IS
     *"/"  CONSIDERED STOPPED.
     *"/" B(I)--ITH COMPONENT OF UNIFORM MAGNETIC FIELD (GAUSS)
     *"/" V--ELECTRIC POTENTIAL VALUE AT THE ORIGIN FOR UNIFORM
     *"/"  POTENTIALS DEFINED AS PHI=-V*(1-Z/ZLIM)
     *"/" STPTIM--TIME WHEN TRAJECTORY ENDS IF IMAGE PLANE NOT HIT")
801   FORMAT(" DTIME--INITIAL INTEGRATION STEP (NS)
     *"/" (NON) UNIFORM POTENTIALS--DISABLE/ENABLE UNIFORM POTENTIAL
     *"/" (NON) UNIFORM B FIELDS--DISABLE/ENABLE UNIFORM B FIELDS
     *"/" TRACK (OFF)--ENABLE/DISABLE POINT BY POINT PRINTOUT OF
     *"/"  PARTICLE TRAJECTORY
     *"//" PARAMETER VALUES MAY BE CHANGED BY SPECIFYING NEW VALUES!
     *"//" STPTIM=5"//" OR BY THROWING A SWITCH!"//" TRACK"//" TRACK
     * OFF"//" ENTER LIST FOR A CURRENT LIST OF PARAMETERS > VALUES
     *"/" ENTER GO TO FINISH")
900   FORMAT(/" YOU MAY CHANGE ANY PARAMETER NOW")
1000  FORMAT(/" BEWARE!! YOUVE CHOSEN A DERIVATIVE STEP SIZE, DELTXI,
     *"/" WHICH MAY CAUSE SAMPLING OF FIELDS OUTSIDE OF THE RADIAL
     *"/" BOUNDARY, RLIM. IF THIS OCCURS, THE PROGRAM WILL EXIT WITH
     *"/" A MESSAGE")
      END
      SUBROUTINE GET(NVALS,VALUES)
C
C THIS ROUTINE GRABS NVALS FREE-FORM VALUES, SEPARATED BY COMMAS.
C
      DIMENSION VALUES(NVALS)
      COMMON/CARD/ICHR,LAST,ICARD(80)
      DATA ICOM/","/
C
1     CALL GETCRD(DUM)
      N=0
2     N=N+1
      IF(ICHR.GT.LAST) GO TO 3
      VALUES(N)=FPBCD(DUM)
      IF(N.EQ.NVALS) RETURN
      IF(ICARD(ICHR).NE.ICOM) GO TO 3
      ICHR=ICHR+1
      GO TO 2
3     PRINT 100
      GO TO 1
100   FORMAT(" ******TOO FEW VALUES SPECIFIED, OR MISSING A COMMA. REENT
     *ER THE LINE")
      END
      SUBROUTINE SCAN(MATCH,NUM)
C
C THIS ROUTINE SCANS THE NEXT FEW CHARACTERS OF THE CURRENT RECORD FOR
C STRINGS WHICH ARE DEFINED IN MATCH. IF THE MATCH IS MADE, THE SERIAL
C NUMBER OF THE STRING IS RETURNED IN NUM, AND THE POINTER BUMPED PAST
C THE STRING.
C
      COMMON/CARD/ICHR,LAST,ICARD(80)
      DIMENSION MATCH(1)
C
      NUM=0
      IEND=1
      NSTRNG=MATCH(1)
      DO 2 M=1,NSTRNG
        IDX=ICHR
        IBGN=IEND+2
        IEND=IBGN+MATCH(IBGN-1)-1
        DO 1 N=IBGN,IEND
          IF(ICARD(IDX).NE.MATCH(N)) GO TO 2
          IDX=IDX+1
1         CONTINUE
        NUM=M
        ICHR = IDX
        RETURN
2       CONTINUE
      RETURN
      END
      SUBROUTINE SCANC(ICHAR)
C
C THIS ROUTINE SCANS THE CURRENT RECORD FOR CHARACTER ICHAR ,STARTING
C AT THE POINTER, AND THEN MOVES THE POINTER JUST PAST IT.
C
      COMMON/CARD/ICHR,LAST,ICARD(80)
C
      IDX=ICHR
      DO 1 I=IDX,LAST
        ICHR=ICHR+1
        IF(ICARD(I).EQ.ICHAR) RETURN
 1      CONTINUE
      ICHR=IDX
      PRINT 100,ICHAR
100   FORMAT(" ******ERROR------LOOKING FOR  "A1)
      RETURN
      END
      SUBROUTINE GETCRD(ANS)
C
C THIS ROUTINE READS A RECORD AND PACKS THE NON-BLANK CHARACTERS INTO
C ARRAY ICARD. IT IGNORES ANY RECORDS STARTING WITH AN ASTRISK (*).
C IT ALSO IGNORES ANY CHARACTERS WHICH PRECEED A DOLLAR SIGN ($) ON
C ANY RECORD.
C
      INTEGER ANS
      COMMON/CARD/ICHR,LAST,ICARD(80)
      DATA IAST,IBLNK,IDOL/"*"," ","$"/
C
      LAST=0
      ICHR=1
C READ IN A CARD
1     READ(5,100) ICARD
      DO 5 I=1,80
      IEND=81-I
      IF(ICARD(IEND).NE.IBLNK) GO TO 7
5     CONTINUE
7     PRINT 200, (ICARD(I),I=1,IEND)
C CHECK FOR A * IN COLUMN 1
      IF(ICARD(1).EQ.IAST) GO TO 1
C REMOVE THE BLANKS FROM THE CARD
      DO 2 I=1,80
        IF(ICARD(I).EQ.IBLNK) GO TO 2
        LAST=LAST+1
        IF(LAST.GT.80) GO TO 6
        ICARD(LAST)=ICARD(I)
2       CONTINUE
C CHECK FOR A $ COMMENT AND SKIP PAST IT.
      DO 3 I=1,LAST
        ICHR=I+1
3       IF(ICARD(I).EQ.IDOL) GO TO 4
      ICHR=1
4     ANS=ICARD(ICHR)
      RETURN
6     PRINT 400
      STOP
100   FORMAT(80A1)
200   FORMAT("     "80A1)
400   FORMAT("  ********TOO MANY CHARACTERS IN CURRENT RECORD-STOP***")
       END
      FUNCTION IBCD(DUMMY)
C
C ROUTINE TO SET UP AN INTEGER FROM THE NEXT FEW CHARACTERS IN THE
C PRESENT RECORD, STARTING AT THE CURRENT POINTER POSITION.
C
      COMMON/CARD/ICHR,LAST,ICARD(80)
      DIMENSION NUMS(10)
      DATA NUMS/"0","1","2","3","4","5","6","7","8","9"/
C
      NUM=0
1     INUM=ICARD(ICHR)
      DO 2 I=1,10
        N=I-1
        IF(INUM.EQ.NUMS(I)) GO TO 3
2       CONTINUE
      GO TO 4
3     NUM=NUM*10+N
      ICHR=ICHR+1
      IF(ICHR.LE.LAST) GO TO 1
4     IBCD=NUM
      RETURN
      END
      FUNCTION FPBCD(DUMMY)
C
C THIS ROUTINE SETS UP A FLOATING POINT QUANTITY FROM THE NEXT SEVERAL
C CHARACTERS IN THE CURRENT RECORD, STARTING AT THE POINTER POSITION.
C
C       CHARACTERS USED........
C            +  OR  -
C            0  THRU  9 (OBVIOUSLY)
C            .   DECIMAL  POINT
C            E   POWER OF 10 FOLLOWS
C            M   TIMES 10**6
C            K   TIMES 10**3
C            U   TIMES 10**-6
C            N   TIMES 10**-9
C            P   TIMES 10**-12
C
      COMMON/CARD/ICHR,LAST,ICARD(80)
      DIMENSION ICHRS(17)
      LOGICAL XPONT,DECML
      DATA ICHRS/"0","1","2","3","4","5","6","7","8","9",
     *     "E","M","K","U","N","P","."/
      DATA IMIN,IPLU/"-","+"/
C
      SIGN=1.
      ISIGN=1
      XPONT=.FALSE.
      DECML=.FALSE.
      XX=0.
      NEXP=0
      DEC=0.
C CHECK 1ST CHAR FOR SIGN
      IT=ICARD(ICHR)
      IF((IT.NE.IMIN).AND.(IT.NE.IPLU)) GO TO 3
      IF(IT.EQ.IMIN) SIGN=-1.
2     ICHR=ICHR+1
      IF(ICHR.GT.LAST) GO TO 16
C LOOP FOR CHARS
3     IT=ICARD(ICHR)
      DO 4 I=1,17
        IDX=I
        IF(IT.EQ.ICHRS(I)) GO TO 5
4       CONTINUE
      GO TO 16
5     GOTO(6,6,6,6,6,6,6,6,6,6,10,11,12,13,14,15,9),IDX
C CHARACTER IS A NUMBER
6     IDX=IDX-1
      IF(XPONT) GO TO 8
      IF(DECML) GO TO 7
C WORKING ON POSITIVE POWERS
      XX=XX*10.+IDX
      GO TO 2
C WORKING ON NEGATIVE POWERS
7     XX=XX+IDX*DEC
      DEC=DEC*.1
      GO TO 2
C WORKING ON THE EXPONENT
8     NEXP=NEXP*10+IDX
      GO TO 2
C CHARACTER IS A DECIMAL POINT
9     IF(DECML.OR.XPONT) GO TO 16
      DEC=.1
      DECML=.TRUE.
      GO TO 2
C CHARACTER IS AN  E
10    XPONT=.TRUE.
C CHECK FOR AN EXPONENT SIGN
      ICHR = ICHR + 1
      IF(ICARD(ICHR).EQ.IPLU) GO TO 2
      IF(ICARD(ICHR).NE.IMIN) GO TO 3
      ISIGN = -1
      GO TO   2
C CHARACTER IS A M
11    NEXP=NEXP+6
      ICHR=ICHR+1
      GO TO 16
C CHARACTER IS A K
12    NEXP=NEXP+3
      ICHR=ICHR+1
      GO TO 16
C CHAR IS A    U
13    NEXP=NEXP-6
      ICHR=ICHR+1
      GO TO 16
C CHAR IS A    N
14    NEXP=NEXP-9
      ICHR = ICHR +1
      GO TO 16
C CHAR IS A    P
15    NEXP=NEXP-12
      ICHR=ICHR+1
C FINISHED........PIECE THE NUMBER TOGETHER
16    FPBCD=SIGN*(XX*10.**(ISIGN*NEXP))
      RETURN
      END
      SUBROUTINE IMAGE
C
C GIVEN THE INITIAL POSITION AND VELOCITY X(I), U(I), OF A CHARGED
C PARTICLE, THIS ROUTINE ACCELERATES THE PARTICLE THROUGH ELECTRO-
C MAGNETIC FIELDS TO IT"S FINAL POSITION ON AN _IMAGE PLANE_.
C T.C. STEPHENS, S.A.I.
C
      DIMENSION UOLD(3),XOLD(3),Y(6),ERR(6)
      LOGICAL TOFINE,BIGHOP,TOCORS,HIT,PAST,HOME
      EXTERNAL DIFFEQ
      COMMON/GUTS/DIAG(5),B(3),V,VUNIF,BUNIF,DELTXI,RLIM,ZLIM,HOPMAX,
     * EP,DIMAG,STPTIM,DTIME
      LOGICAL VUNIF,DIAG,BUNIF
      COMMON/IMAG/ X(3),U(3),T,NSTEP,XERR,YERR,ZERR
      EQUIVALENCE (Y,X)
C
C INITIALIZE VARIABLES.
C
      XERR=0.
      YERR=0.
      ZERR=0.
      ERR(1)=0.
      ERR(2)=0.
      ERR(3)=0.
      NSTEP=0
      TNS=0.
      T=0.
      TOCORS=.FALSE.
      TOFINE=.FALSE.
      BIGHOP=.FALSE.
      HIT=.FALSE.
      PAST=.FALSE.
      HOME=.FALSE.
      DT=DTIME*1.E-9
      IF(DIAG(1)) PRINT 100,T,X,U,DT
C
C INTEGRATE EQUATIONS OF MOTION] MOVE PARTICLE THROUGH FIELDS.
C
1     TOLD=T
      XERR=XERR+ERR(1)*1.E4
      YERR=YERR+ERR(2)*1.E4
      ZERR=ZERR+ERR(3)*1.E4
      DO 2 I=1,3
        UOLD(I)=U(I)
2       XOLD(I)=X(I)
      IF(TNS.GE.STPTIM) RETURN
3     IF(DIAG(5)) PAUSE
      NSTEP=NSTEP+1
      IF(NSTEP.GT.1000) PAUSE"PARTICLE HAS FAILED TO HIT IMAGE PLANE
     *IN 1000 STEPS"
      CALL KUTMER(6,DT,T,Y,DIFFEQ,ERR)
C
C CHECK THAT PARTICLE IS NOT MOVING AWAY FROM THE IMAGE PLANE.
C
      TNS=T*1.E9
      IF(DIAG(1)) PRINT 200,TNS,X,U,DT
      IF(X(3).LT.-1.)STOP"IMAGE!PARTICLE MOVING AWAY FROM IMAGE PLANE"
C
C CHECK THAT THE TIMESTEP IS VALID BY MAKING SURE THAT! 1) THE ERROR IN
C THE COORDINATE INTEGRATIONS IS SMALL 2) THE PARTICLE IS _SEEING_ THE
C GEOMETRY (I.E. NOT MOVING TOO FAR IN ONE STEP) 3) THE TIME STEP MIGHT
C BE INCREASED FOR SPEED 4) THE IMAGE PLANE IS NOT PASSED BEYOND.
C
      ERRMAX=AMAX1(ERR(1),ERR(2),ERR(3))*1.E4
      TOCORS=(ERRMAX.GT.EP)
      TOFINE=(ERRMAX.LT.EP/16.)
      VEL=SQRT(U(1)*U(1)+U(2)*U(2)+U(3)*U(3))
      BIGHOP=((DT*VEL).GT.HOPMAX)
      PAST=(X(3).GT.ZLIM)
      HIT=(ABS(X(3)-ZLIM).LE.DIMAG)
C
      IF(HIT.AND..NOT.(TOCORS.OR.BIGHOP)) RETURN
      T1=DT
      T2=DT
      T3=DT
      IF(TOCORS) T1=DT/2.
      IF(BIGHOP) T2=HOPMAX/VEL/2.
      IF(PAST) T3=(ZLIM-XOLD(3))/U(3)
      IF(PAST) HOME=.TRUE.
      DT=AMIN1(T1,T2,T3)
      IF(.NOT.TOFINE.OR.BIGHOP.OR.PAST.OR.HOME) GO TO 4
      DT=DT*1.1
      IF(DIAG(1)) PRINT 400,DT
      GO TO 1
4     IF(TOCORS.OR.BIGHOP.OR.PAST) GO TO 6
      IF(HOME) DT=(ZLIM-X(3))/U(3)
      GO TO 1
C
C TIME STEP WAS DEFECTIVE] RETURN PARTICLE TO PREVIOUS STATE > TRY AGAIN
C
6     T=TOLD
      NSTEP=NSTEP-1
      IF(.NOT.DIAG(1)) GOTO 7
      IF(TOCORS) PRINT 300, ERRMAX,T1
      IF(BIGHOP) PRINT 500, T2
      IF(PAST) PRINT 600, T3
7     DO 8 I=1,3
        U(I)=UOLD(I)
8       X(I)=XOLD(I)
      GO TO 3
C
100   FORMAT(" TIME(NS)",7X,"X",8X,"Y",8X,"Z",8X,"VX",7X,"VY",7X,"VZ"
     *,8X,"DT"/F8.2,2X,3F9.4,2X,3E9.2,2X,E9.3)
200   FORMAT(F8.2,2X,3F9.4,2X,3E9.2,2X,E9.3)
300   FORMAT(8X,"ONE COORD HAD ERROR OF "F7.2," MICRONS...CUT DT TO"
     *E9.3)
400   FORMAT(8X,"ALL COORD ERRORS LESS THAN EP/16 MICRONS...BUMP DT UP
     * TO"E9.3)
500   FORMAT(8X,"PARTICLE MOVED TOO FAR LAST STEP] CUT DT TO"E9.3)
600   FORMAT(8X,"PARTICLE CROSSED IMAGE PLANE] BACKUP] DT="E9.3)
      END
      SUBROUTINE DIFFEQ(T,Y,DYDT)
C
C THIS CODE PROVIDES THE FORM OF THE DIFFERENTIAL EQUATIONS FOR A
C CHARGED PARTICLE IN AN ELECTRO-MAGNETIC FIELD. IT IS TAYLORED
C ESPECIALLY TO THE KUTTA-MERSON INTEGRATOR, KUTMER.
C T.C. STEPHENS S.A.I.
C
      DIMENSION Y(6),DYDT(6),U(3),X(3),DPDX(3)
      COMMON/GUTS/DIAG(5),B(3),V,VUNIF,BUNIF
      LOGICAL VUNIF,DIAG,BUNIF
      DATA C,E,EM/2.997925E10,-4.803250E-10,9.109558E-28/
C
      DO 1 I=1,3
        X(I)=Y(I)
        U(I)=Y(I+3)
1       DYDT(I)=U(I)
C
      DO 3 I=1,3
        IX=I
3       CALL DERIV(IX,X,PHI,DPDX(I))
      IF(.NOT.BUNIF) CALL NUBFLD(1,X,B)
      VEL=SQRT(U(1)*U(1)+U(2)*U(2)+U(3)*U(3))
      IF(VEL.GT.C) STOP"DIFFEQ! TIME STEP GUESS IS TOO LARGE] CUT BACK"
      GAM=1./SQRT(1.-(VEL/C)**2)
      UDDELP=U(1)*DPDX(1)+U(2)*DPDX(2)+U(3)*DPDX(3)
C
      DYDT(4)=E/(GAM*EM)*(-DPDX(1)+U(2)*B(3)/C-U(3)*B(2)/C+U(1)*UDDELP
     * /C**2)
      DYDT(5)=E/(GAM*EM)*(-DPDX(2)+U(3)*B(1)/C-U(1)*B(3)/C+U(2)*UDDELP
     * /C**2)
      DYDT(6)=E/(GAM*EM)*(-DPDX(3)+U(1)*B(2)/C-U(2)*B(1)/C+U(3)*UDDELP
     * /C**2)
C
      IF(DIAG(3)) PRINT 100,X,U,DYDT
100   FORMAT(5X,"DIFFEQ!X="3F8.4," U="3E9.3/,5X,"DYDT="6E9.3)
      RETURN
      END
      SUBROUTINE DERIV(IX,X,PHI,DPDXI)
C
C THIS ROUTINE TAKES THE DERIVATIVE OF THE SCALAR POTENTIAL,
C PHI, AT POINT X(I) (TO THIRD ORDER), W.R.T. THE IX"TH CO-ORD.
C DPHI/DX(IX) IS PUT IN DPDXI.
C T.C. STEPHENS S.A.I.
C
      DIMENSION XX(3),VV(4),X(3),COEF(4,4),XYZ(3)
      COMMON/GUTS/DIAG(5),B(3),V,VUNIF,BUNIF,DELTXI,RLIM,ZLIM
      LOGICAL VUNIF,DIAG,BUNIF
      DATA XYZ/"X","Y","Z"/
      DATA COEF/-11.,-2.,1.,-2.,18.,-3.,-6.,9.,-9.,6.,3.,-18.,
     * 2.,-1.,2.,11./
C
C GIVEN PROBLEM BOUNDARIES > CURRENT POINT OF INTEREST, SELECT WHICH OF
C THE FOUR CONSECUTIVE POINTS (USED TO DETERMINE 3RD ORDER DERIVATIVES)
C WILL BE THE POINT AT WHICH THE DERIVATIVE IS EVALUATED (I.E. CHOOSE
C DIRECTION OF SAMPLING).
C
      IPNT=2
      IF(IX.NE.3) GO TO 6
C
C Z DIRECTION DECISION...
C
      IF(X(3).LT.DELTXI) IPNT=1
      IF(X(3).GT.(ZLIM-2.*DELTXI)) IPNT=4
      GO TO 7
C
C X OR Y DIRECTION DECISION...
C
6     IF((IX.EQ.1).AND.(X(1).GT.0.)) IPNT=3
      IF((IX.EQ.2).AND.(X(2).GT.0.)) IPNT=3
C
C GET SCALAR POTENTIAL VALUES AT FOUR POINTS (STORE IN VV).
C
7     DO 1 I=1,3
1       XX(I)=X(I)
      DO 3 J=1,4
        XX(IX)=(J-IPNT)*DELTXI+X(IX)
        IF(SQRT(XX(1)**2+XX(2)**2).GT.RLIM)STOP"DERIV!SAMPLING OUTSIDE
     *  OF PROBLEM"
        IF(VUNIF) PHI=-V*(1.-(XX(3)/ZLIM))
        IF(.NOT.VUNIF) CALL NUPOT(1,XX,PHI)
        VV(J)=PHI
3       CONTINUE
C
C USE THE 4 POINT LEGENDRE DERIVATIVE FORMULA FROM ABRAMOWITZ > SETGUN,
C PP 914 WHERE IPNT IS  THE POINT OF THE DERIVATIVE.
C
      V1=VV(1)
      V2=VV(2)
      V3=VV(3)
      V4=VV(4)
      DPDXI=(COEF(IPNT,1)*V1+COEF(IPNT,2)*V2+COEF(IPNT,3)*V3
     * +COEF(IPNT,4)*V4)/(6.*DELTXI)
      PHI=VV(IPNT)
C
      IF(DIAG(4)) PRINT 100,XYZ(IX),X,IPNT,DELTXI,DPDXI,PHI
100   FORMAT(6X,"DERIV! W.R.T. "A1," AT X,Y,Z="3F8.4," IPNT,DELTXI="I2,
     *F8.4/6X,"DPDXI,PHI= "2E9.3)
      RETURN
      END
      SUBROUTINE KUTMER(N,H,X,Y,FXY,E)
C
C THIS ROUTINE INTEGRATES SETS OF FIRST ORDER ORDINARY DIFFERENTIAL
C EQUATIONS BY THE KUTTA-MERSON METHOD (SEE FOX, _NUMERICAL SOLUTION
C OF ORDINARY > P.D.E._, PP24).
C    N= # OF EQUATIONS
C    H= STEPSIZE
C    X= INDEPENDENT VARIABLE VALUE
C    Y= ARRAY OF DEPENDENT VARIABLE VALUES
C    FXY= USER DEFINED SUBROUTINE WHICH DEFINES THE DERIVATIVES
C    E= ARRAY OF ERROR ESTIMATES IN THE RESPECTIVE DEPENDENT VARIABLES
C
      DIMENSION F0(10),F2(10),F3(10),YZERO(10),Y5(10),Y(10),DY(10),E(10)
      EQUIVALENCE(Y5(1),F2(1))
      IF(N.GT.10) STOP"KUTMER!TOO MANY DEPENDENT VARIABLES"
      DO 1 I=1,N
1     YZERO(I)=Y(I)
      CALL FXY (X,Y,DY)
      DO 2 I=1,N
      F0(I) = H*DY(I)
2     Y(I)=YZERO(I)+.333333333333333*F0(I)
      X=X+.333333333333333*H
      CALL FXY (X,Y,DY)
      DO 3 I=1,N
3     Y(I)=YZERO(I)+.166666666666667*F0(I)+.166666666666667*H*DY(I)
      CALL FXY (X,Y,DY)
      DO 4 I=1,N
      F2(I)=H*DY(I)
4     Y(I)=YZERO(I)+.125*F0(I)+.375*F2(I)
      X=X+.166666666666667*H
      CALL FXY (X,Y,DY)
      DO 5 I=1,N
      F3(I)=H*DY(I)
5     Y(I)=YZERO(I)+.5*F0(I)-1.5*F2(I)+2.*F3(I)
      X=X+.5*H
      CALL FXY (X,Y,DY)
      DO 6 I=1,N
6     Y5(I)=YZERO(I)+.166666666666667*F0(I)+.666666666666667*F3(I)
     *+.166666666666667*H*DY(I)
7     DO 8 I=1,N
      E(I)=.2*ABS(Y5(I)-Y(I))
      Y(I)=Y5(I)
8     CONTINUE
      RETURN
      END
      SUBROUTINE SOURCE(MODE,Q,V1,V2,V3)
C
C THIS IS A SAMPLE MONTE CARLO SOURCE OF ELECTRONS WHICH PRODUCES
C PARTICLES WITH AN ENERGY DISTRIBUTION F(E)=6.*(EM*E-E**2)/EM**3
C WHERE EM IS THE MAXIMUM ENERGY OF THE SPECTRUM (EV). THE ANGLE
C DISTRIBUTIONS ARE UNIFORM AZIMUTHALLY AND COS**2 ABOUT THE Z AXIS.
C WHEN CALLED WITH MODE=1, THIS ROUTINE RETURNS WITH A PARTICLE WEIGHT,
C Q, AND VELOCITY COMPONENTS, V1,V2,V3. T.N. DELMER
C
      F(E)=6./EM**3*(EM*E-E**2)
C
      IF(MODE.NE.0) GO TO 1
      PRINT 100
100   FORMAT(" SUBROUTINE SOURCE(MODE,Q,V1,V2,V3) HAS JUST BEEN CALLED
     *"/" WITH MODE=0 THIS IS THE INITIALIZATON PHASE FOR THE USER-
     *"/" DEFINED MONTE CARLO SOURCE. IN THIS CASE, THE ENERGY
     *"/" DISTRIBUTION FOR ELECTRONS IS!  F(E)=6.*(EM*E-E**2)/EM**3
     *"/" ANGLE DISTRIBUTONS ARE UNIFORM AZIMUTHALLY] COS**2 IN SOLID
     *"/" ANGLE ABOUT THE Z AXIS. NOW ENTER EM, THE MAXIMUM ENERGY OF
     *"/" THE SPECTRUM, IN EV.")
      CALL GET(1,EM)
      RETURN
1     X=0.
      XP=RANDOM(X)
      XT=RANDOM(X)
      XE=RANDOM(X)
      PH=2.*3.14159*XP
      TH=ATAN2(SQRT(1.-XT**2),XT)
C THE ENERGY IS SELECTED FLAT FROM ZERO TO EM
      E=EM*XE
      Q=F(E)*EM*(COS(TH))**2*3.
C THIS CHOICE OF THE WEIGHTING FUNCTION COMPENSATES FOR THE FLAT ENERGY
C SELECTION TO REPRODUCE THE  SPECIFIED SPECTRUM. THIS ANGLE WEIGHTING
C GENERATES A COS**2 DISTRIBUTION. IF F IS NORMALIZED, THE AVERAGE Q IS
C UNITY
      V=SQRT(2.*E*1.6E-12/9.11E-28)
      V1=V*SIN(TH)*COS(PH)
      V2=V*SIN(TH)*SIN(PH)
      V3=V*COS(TH)
      RETURN
      END
      FUNCTION RANDOM(DUM)
C
C DUMMY RANDOM NUMBER GENERATOR (UNIFORM)
C
      RANDOM=0.
      PRINT 100
100   FORMAT(" YOU MUST SUPPLY A UNIFORM RANDOM NUMBER GENERATOR FOR THI
     *S SOURCE OF "/" ELECTRONS] THEY ARE USUALLY MACHINE DEPENDENT")
      STOP"RANDOM!NEED A RANDOM NUMBER GENERATOR"
      END
      SUBROUTINE NUBFLD(MODE,X,B)
C
C THIS IS A SAMPLE USER-DEFINED NON-UNIFORM B FIELD. IT SIMULATES A
C HELMHOLTZ PAIR CENTERED ON THE Z AXIS.
C
      DIMENSION X(3),B(3),XINP(4)
      IF(MODE.NE.0) GO TO 1
      PRINT 100
100   FORMAT(" SUBROUTINE NUBFLD(MODE,X,B) HAS JUST BEEN CALLED WITH
     *"/" MODE=0. THIS IS THE INITIALIZATION PHASE FOR THE USER-DEFINED
     *"/" NON-UNIFORM B FIELD, WHICH IN THIS SAMPLE CASE SIMULATES THE
     *"/" MAGNETIC FIELD NEAR A HELMHOLTZ PAIR OF COILS CENTERED ON THE
     *"/" TUBE AXIS. YOU MUST NOW SPECIFY THE DESIRED FIELD STRENGTH
     *"/" (GAUSS) ON THE AXIS, THE COIL RADIUS (CM), THE TUBE RADIUS
     *"/" (RLIM IN CM), AND THE Z POSITION OF THE CENTER OF THE COILS")
      CALL GET(4,XINP)
      BZ0=XINP(1)
      A=XINP(2)
      RLIM=XINP(3)
      ZCNTR=XINP(4)
      D=A/2.
      BCOEF=BZ0*A**3*5.*SQRT(5.)/16.
      ORDER=(RLIM/A)**4
      PRINT 200,ORDER
200   FORMAT(" FIELDS PRODUCED ARE ACCURATE TO ORDER (R/A)**4, WHERE
     *"/" R= RADIAL POSITION OF ELECTRON] A= COIL RADIUS. AT THE TUBE
     *"/" EDGE (RLIM), ACCURACY IS OF ORDER "G10.4)
      RETURN
1     Z=X(3)-ZCNTR
      R1=SQRT((Z+D)**2+A**2)
      R2=SQRT((Z-D)**2+A**2)
      R=SQRT(X(1)**2+X(2)**2)
      B(3)=BCOEF*(1./R1**3+1./R2**3+.75*R**2*((A**2-4*(Z+D)**2)/R1**7+
     * (A**2-4*(Z-D)**2)/R2**7))
      TERM1=((Z+D)/R1**5+(Z-D)/R2**5)
      TERM2=((3*A**2*(Z+D)-4*(Z+D)**3)/R1**9+(3*A**2*(Z-D)-4*(Z-D)**3)/
     *  R2**9)
      BR=BCOEF*(1.5*R*TERM1+15./16.*R**3*TERM2)
      THE=ATAN2(X(2),X(1))
      B(1)=BR*COS(THE)
      B(2)=BR*SIN(THE)
      RETURN
      END
      SUBROUTINE NUPOT(MODE,X,PHI)
C
C THIS IS A SAMPLE USER-DEFINED NON-UNIFORM POTENTIAL. IT PRODUCES AN
C INVERSE SQUARE FORCE FIELD ABOUT THE Z AXIS (F=-A/R**2 ] V=-A/R)
C
      DIMENSION X(3)
      IF(MODE.NE.0) GO TO 1
      PRINT 100
100   FORMAT(" SUBROUTINE NUPOT(MODE,X,PHI) HAS JUST BEEN CALLED WITH
     *"/" MODE=0 THIS IS THE INITIALIZATION PHASE FOR THE USER-DEFINED
     *"/" NON-UNIFORM POTENTIAL, WHICH IN THIS SAMPLE CASE PRODUCES AN
     *"/" INVERSE SQUARE FIELD ABOUT THE Z AXIS
     *"/"  F=-A/R**2 ] V=-A/R
     *"/" NOW, YOU MUST ENTER THE FIELD STRENGTH, A,(STATVOLT-CM)!")
      CALL GET(1,A)
      RETURN
1     R=SQRT(X(1)**2+X(2)**2)
      IF(R.EQ.0.)STOP"NUPOT!ELECTRON IS ON Z AXIS! INFINITE POTENTIAL"
      PHI=A/R
      RETURN
      END
